import numpy as np

import sympy
from fractions import Fraction
np.set_printoptions(formatter={'all':lambda x: str(Fraction(x).limit_denominator())})#设置np.array用分数表示
A = np.array([[1,1,2,4], [-1,-1,-1,-2], [3,3,2,4]])
b = np.array([4,1,4]).reshape(-1,1)#创建列向量

def leastNorm(A,b):
    """
    Ax=b
    求：（1）广义逆（也叫“伪逆”）；（2）极小范数最小二乘解;(3)相容判断(满足则相容)：r(A)==r(A,b)
    (4)最小二乘解=极小范数最小二乘解+(I-A+A)t)
    """
    #（1）广义逆（也叫“伪逆”）
    A_reve_global=np.linalg.pinv(A)#广义逆（也叫“伪逆”）-非方阵求逆
    print("广义逆（也叫“伪逆”）A+:\n",A_reve_global)
    #Ax=b的极小范数最小二乘解为A+b

    #（2）极小范数最小二乘解
    least_norm=A_reve_global.dot(b)
    print("Ax=b的极小范数最小二乘解为A+b:\n",least_norm)

    #(3)相容判断(满足则相容)：r(A)==r(A,b)
    Ab=np.hstack((A,b))#方程的增广矩阵-矩阵横向连接
    # print("Ab=",Ab)
    rA=np.linalg.matrix_rank(A)  # 矩阵的秩
    rAb=np.linalg.matrix_rank(Ab)#增广矩阵的秩
    # print("rA=",rA,"rAb=",rAb)
    flag = (rA == rAb)
    print("Ax=b相容判断:",flag )

    #(4)最小二乘解
    ApA=A_reve_global.dot(A)#A+A
    # print(ApA)
    IApA=np.eye(ApA.shape[0])#I单位矩阵
    # print(IApA)
    x, y, t = sympy.symbols('x y t')#符号
    least_square = least_norm + (IApA - ApA)*t
    print("Ax=b的极小范数最小二乘解为A+b:\n", least_square)

leastNorm(A,b)